#include<iostream>
#include<cmath>
#include<vector>
#include "../programming1/EquationSolver.h"
#include "../project1/LUFactorization.h"
using namespace std;

class CubicSpline {
private:
	vector <double> x, f, m, M;	
public:
    CubicSpline(const vector<double>& x, const vector<double>& f, const string& mode = "natural", const double& m_0 = 0, const double& m_n = 0) : x(x), f(f) {
		int n = x.size() - 1;
		m.resize(n+1);
		vector <double> lamda(n+1), mu(n+1);	
		for (int i = 1; i <= n-1; ++ i) {
			lamda[i] = (x[i] - x[i-1]) / (x[i+1] - x[i-1]);
			mu[i] = (x[i+1] - x[i]) / (x[i+1] - x[i-1]);
		}
    
        vector <double> f1(n+1), f2(n+2);	// 一阶差商, 二阶差商
        vector <double> d(n+1), u(n), l(n), b(n+1);

        if (mode == "natural" && (m_0 != 0 || m_n != 0)) throw "The boundary condition is wrong!";
        if (mode == "complete" || mode == "natural"){
            for (int i = 1; i <= n; i++)
            {
                f1[i] = (f[i] - f[i-1]) / (x[i] - x[i-1]);
            }
            f2[1] = (f1[1] - m_0) / (x[1] - x[0]);
            for (int i = 2; i <= n; ++ i){
                f2[i] = (f1[i] - f1[i-1]) / (x[i] - x[i-2]);
            }
            f2[n+1] = (m_n - f1[n]) / (x[n] - x[n-1]);

            for (int i = 1; i <= n-1; ++ i) {
                    d[i] = 2; //对角
                    u[i-1] = mu[i]; //上面的对应/mu
                    l[i] = lamda[i]; //下面的对应/lamda
                    b[i] = 6 * f2[i+1]; 
                }
            d[0] = 2, u[0] = 1, b[0] = 6 * f2[1];
            d[n] = 2, l[n-1] = 1, b[n] = 6 * f2[n+1];
            // 调用LU分解求解三对角方程组
            M=LU(d,u,l,b);
            m[0] = m_0, m[n] = m_n;
            for (int i = 1; i <= n-1; ++ i){
                m[i] = f1[i+1] - (2 * M[i] + M[i+1]) * (x[i+1] - x[i]) / 6;
            }		
        }
        else if(mode == "ssd"){
            for (int i = 1; i <= n; ++ i)
                    f1[i] = (f[i] - f[i-1]) / (x[i] - x[i-1]);
                for (int i = 2; i <= n; ++ i)
                    f2[i] = (f1[i] - f1[i-1]) / (x[i] - x[i-2]);
                for (int i = 1; i <= n-1; ++ i) {
                    d[i] = 2;
                    u[i-1] = mu[i];
                    l[i] = lamda[i];
                    b[i] = 6 * f2[i+1];
                }
                d[0] = 1, b[0] = m_0;
                d[n] = 1, b[n] = m_n;
                M = LU(d, u, l, b);
                for (int i = 0; i <= n-1; ++ i)
                    m[i] = f1[i+1] - (2 * M[i] + M[i+1]) * (x[i+1] - x[i]) / 6;
                m[n] = f1[n] - (2 * M[n] + M[n-1]) * (x[n-1] - x[n]) / 6;
        }
    }       
    double getValue(const double& _x) const {
		int i = upper_bound(x.begin(), x.end(), _x) - x.begin() - 1;
		double h = _x - x[i], c0 = f[i], c1 = m[i], c2 = M[i] / 2, c3 = (M[i+1] - M[i]) / (x[i+1] - x[i]) / 6;
		return c0 + h * (c1 + h * (c2 + c3 * h));

    }
};

CubicSpline  CubicSplineInterpolation(Function & f,  const double& l, const double& r, const int& n, const string& mode = "natural") {
	vector <double> x(n+1), y(n+1);
	for(int i = 0; i <= n; ++ i) x[i] = l + (r - l) / n * i, y[i] = f(x[i]);
	double m_0, m_n;
	if (mode == "natural") m_0 = m_n = 0;
	else if (mode == "complete") m_0 = f.diff(l), m_n = f.diff(r);
	else if (mode == "ssd") m_0 = f.diff2(l), m_n = f.diff2(r);
	return CubicSpline(x, y, mode, m_0, m_n);
}